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Interpolating scaling functions give a faithful representation of a localized charge distribution by 
its values on a grid. For such charge distributions, using a Fast Fourier method, we obtain highly 
accurate electrostatic potentials for free boundary conditions at the cost of 0(N logN) operations, 
where N is the number of grid points. Thus, with our approach, free boundary conditions are 
treated as efficiently as the periodic conditions via plane wave methods. 
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I. INTRODUCTION 

Solving Poisson's equation 
V 2 V = 



(1) 



to find the electrostatic potential V arising from a charge 
distribution p is a basic problem that can be found in 
nearly any field of physics and chemistry. It is therefore 
essential to have efficient solution methods for it. 

A large variety of methods has been developed for sys- 
tems of point particles interacting by electrostatic forces. 
Formally this problem can be considered as the solution 
of Poisson's equation where the charge distribution is a 
sum of delta functions. The classical method for periodic 
boundary conditions is the Ewald method (1). For large 
systems and free boundary conditions the Fast Multipole 
Method (2) is a powerful method due to its linear scaling 
with respect to the number of particles. 

The FMM method has been generalized for continous 
systems where the charge density can be represented as 
a sum of Gaussians multiplied by a polynomial (3). This 
generalization exhibits linear scaling with respect to the 
volume but not with respect to the number of Gaus- 
sians at constant volume. It works well in the context of 
quantum chemistry calculations with medium size Gaus- 
sian basis sets where the charge density is naturally ob- 
tained in the required form but it is not a general purpose 
method. 

For periodic boundary conditions and smoothly vary- 
ing charge densities plane wave methods are simple and 
powerful because the Laplacian is a diagonal matrix in 
a plane wave representation. Given a Fourier represen- 
tation of the charge density one can therefore obtain the 
potential simply by dividing each Fourier coefficient by 
Ikl 2 , where k is the wavevector of the Fourier coefficient. 



If the charge density is originally given in real space a first 
Fast Fourier Transformation (FFT) is required to obtain 
its Fourier coefficients, and a second FFT is required to 
obtain the potential in real space from its Fourier coeffi- 
cients. The overall scaling is therefore of order NlogN 
where N is the number of grid points. 

Because of the simplicity of these plane wave methods 
various attempts have been made to generalize them to 
free boundary conditions. The most rudimentary method 
is just to take a very large periodic volume and to hope 
that the amplitude of the potential is nearly zero on the 
surface of the volume. Due to the long range of elec- 
trostatic forces this condition is however not fulfilled for 
volume sizes that are affordable with plane waves. Such 
as scheme is only possible if adaptive periodic wavelets 
are used (6). In addition periodic boundary conditions do 
not permit to treat systems with monopoles and dipoles 
because for such systems no well defined solution exists 
under periodic boundary conditions. The first method 
to attack the problem in a systematic way was by Hock- 
ney (4). He proposed a Fourier approximation to the 
kernel 



K(r) = - 



(2) 
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of the Poisson equation. The method was intended for 
applications in plasma physics where no high accuracy is 
required. In other application such as electronic structure 
calculations high accuracy is required and the method is 
not optimal. For a spherical geometry the Fourier coeffi- 
cients of the i kernel can be calculated analytically. This 
is the basis of the simple and powerful method by Fiisti- 
Molnar and Pulay (11). Its obvious restriction is that it 
is efficient only for spherical geometries. Another rather 
complicated method was proposed by Martyna and Tuck- 
erman. This method gives high accuracy only in the cen- 
ter of the computational volume. It requires therefore ar- 
tificially large simulation boxes which is numerically very 
expensive. All the above discussed methods use FFT's 
at some point and have therfore a 0(N log N) scaling. 



2 



In this paper we will describe a new Poisson solver for 
free boundary conditions on an uniform mesh. Contrary 
to Poisson solvers based on plane wave functions, our 
method is using interpolating scaling functions to repre- 
sent the charge density. It is therefore from the begin- 
ning free of long range interactions between supcrcclls, 
that falsify results if plane waves are used to describe 
non-periodic systems. Due to the convolutions we have 
to evaluate our method has a N log(iV) scaling instead 
of the ideal linear scaling, Due to its small prefactor the 
method is however most efficient when dealing with lo- 
calized densities such as can be found for example in the 
context of ab initio pseudo-potential electronic structure 
calculations using finite differences (8) finite elements (9) 
or plane waves for non-periodic systems. 



II. INTERPOLATING SCALING FUNCTIONS 

Scaling functions arise in wavelet theory (5). A scaling 
function basis set can be obtained from all the transla- 
tions by a certain grid spacing h of the mother wavelet 
centered at the origin. What distinguished scaling func- 
tions from other localized basis functions like finite ele- 
ments is the refinement relation. The refinement relation 
establishes a relation between a scaling function <p(x — i) 
and the same scaling functions compressed by a factor of 
two, or, equivalently, between the scaling functions on a 
grid with grid spacing h and another one with spacing 
h/2. For the scaling function centered a the origin, it 
reads 
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FIG. 1 Plots of interpolating scaling functions and wavelets 
of 14-th and 100-th order. 



4>{x)= h j( j>(2x-j) 



(3) 



where the hj 's are the elements of a filter that character- 
izes the wavelet family. An interpolating wavelet has in 
addition the property that it is equal to one at the ori- 
gin and zero at all other integer points. Because of this 
property it is very simple to find the scaling function 
expansion coefficients of any function. The coefficients 
are just the values of the function to be expanded on 
the grid, m-th order interpolating scaling functions are 
generated by m — 1-th order recursive interpolation (7). 
Fig. 1 shows an 14-th order and 100-th order interpolat- 
ing scaling function. Three-dimensional scaling functions 
can be obtained as the product of their one-dimensional 
counterparts 

$ ii,i2,i 3 ( r ) = <f>ii(x) <t>i 2 {y) <f>i 3 (z) 

= <j>{x - h) 4>{y - i 2 ) 4>{z - i 3 ) (4) 

where r = (x,y, z). The points ii,fl2,*3 are the nodes 
of an uniform 3-dimensional mesh, with i p = 1, • • • , n p , 
P= 1,2,3. 

Continuous charge distributions are represented in nu- 
merical work typically by their values Pij,k on a grid. It 
follows from the above described properties of interpo- 
lating scaling functions that the corresponding continous 
charge density is given by 

P( r ) = E Pn,i2,i3<t>( x - *i) <t>(y - *2) 4>{z - h) (5) 



The mapping of Eq. 5 between the discretized and conti- 
nous charge distribution ensures that the first m discrete 
and continous moments are identical for a m-th order 
interpolating wavelet family, i.e. 



i,j,k 



':,3,k = J < 



dr x ei y t2 z l3 p(r) 



(6) 



if ^1,^2,^3 < Tn. The proof of this relation is given in 
the Appendix. Since the various multipoles of the charge 
distribution determine the major features of the poten- 
tial the above equalities tell us that a scaling function 
representation gives the most faithful mapping between 
a continuous and discretized charge distribution for elec- 
trostatic problems. 



III. POISSON'S EQUATION IN A BASIS SET OF 
INTERPOLATING SCALING FUNCTIONS 

As is well known the following integral equation gives 
the potential for free boundary conditions 



V(r) 



I 



dr'- 



1 



(7) 



We are interested in the values of the potential on the 
same grid that was used for the charge density. Denoting 
the potential on the grid point rj 1 j 2 j 3 = (xj 1 , yj 2 , Zj 3 ) by 
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V h, 32,33 = V ( T h, 32,33) we have 



(8) 



dr 



,<i) ll {x') <p t2 (y') <f> l3 (z') 



I\,l2,l3 

The above integral defines the discrete kernel 
dr> n (x')<My') M*'h — 



(9) 



Since the problem is invariant under combined transla- 
tions of both the source point (ii, 12,13) and the obser- 
vation point (ii, 72,73) the kernel depends only on the 
difference of the indices 

K{i\,j\;i2,h\is,h) = K(h -31,12 -32,13 - h) (10) 

and the potential Vj 1 j 2 j 3 can be obtained from the 
charge density Pi x ,i 2 ,i 3 by the following 3-dimensional 
convolution: 



K(i 1 -j 1 ,i 2 -j 2 ,i3-j3)pi 1 ,i 2 ,i 3 . (11) 



Once the Kernel is available in Fourier space, this con- 
volution can be evaluated with two FFTs at a cost of 
0(N log N) operations where N = n\ n 2 ^3 is the num- 
ber of 3-dimensional grid points. Since all the quantities 
in the above equation are real, rcal-to-complcx FFT's can 
be used to reduce the number of operations compared to 
the case where one would use ordinary complex-complex 
FFT's. Obtaining the Kernel in Fourier space from the 
Kernel K(j\, 32,33) m real space requires another FFT. 

It remains now to calculate the values of all the 
elements of the kernel K(ki,k 2 ,k 3 ). Solving a 3- 
dimcnsional integral for each element would be much too 
costly and we use therefore a separable approximation of 
1/r in terms of Gaussians (12; 13), 



(12) 



In this way all the complicated 3-dimensional integrals 
become products of simple 1-dimensional integrals. Us- 
ing 89 Gaussian functions with the coefficients u k and pk 
suitably chosen, we can approximate - with an error less 
than 10~ 8 in the interval [10~ 9 ,1]. If we are interested 
in a wider range, e.g. a variable R going from zero to L, 



we can use r = ^: 

L 
R 
1 

R 



k 



Pk f>2 



L 



uj k e 



-P k R 2 



Pk = ^ 



Pk_ 

L 2 



With this approximation, we have that 



89 

E 



io h K h (p k )K h {p k )K h (p k ) 



(13) 
(14) 
(15) 

(16) 



where 

Kj(p k ) = J ^(x)e-P^ 2 dx (17) 
= J <p (x)e- pk( > x -rf dx . (18) 

So we only need to evaluate 89 x max ({ni, n 2 , n 3 }) in- 
tegrals of the type 



KM 



(x)e- p{x -rt 2 dx , 



(19) 



for some value of p chosen between 3 • 10~ 5 and 3 • 10 16 . 

The accuracy in calculating the integrals can be further 
improved by using the refinement relation for interpolat- 
ing scaling functions 3. 

From 19, we can evaluate Ki(Ap) as: 

Ki(Ap) = J <p(x)e- 4p{x - l)2 dx (20) 
= i J V {x/2)e- p{x - 2l ^dx (21) 

= \Y. h ij ^)e- p{x ~ 2i)2 dx (22) 
3 

= \ Y. h ^-M ■ (23) 



The best accuracy in evaluating numerically the integral 
is attained for p < 1. For a fixed value of p given by 
Eq. 12, the relation 23 is iterated n = [log 4 (p)] times 
starting with p = So the numerical calculation of 
the integrals Ki(j>) is performed as follows: for each p, we 
compute the number n of required recursions levels and 
calculate the integral Ki(p ). The value of n is chosen 
such that po ~ 1 so we have a gaussian functions not 
too sharp. The evaluation of the interpolating scaling 
functions is fast on a uniform grid of points so we perform 
a simple summation over all the grid points. In Figure 2, 
we show that 1024 points are enough to obtain machine 
precision. Note that the values of K (p) vary over many 
orders of magnitude as shown in In Figure 2. 



IV. NUMERICAL RESULTS AND COMPARISON WITH 
OTHER METHODS 

We have compared our method with the plane wave 
methods by Hockncy (4) and Martyna and Tucker- 
man (10) as implemented in the CPMD electronic struc- 
ture program (15). As expected Hockney's method does 
not allow to attain high accuracy. The method by Mar- 
tyna and Tuckerman has a rapid exponential convergence 
rate which is characteristic for plane wave methods. Our 
new method has an algebraic convergence rate of h m with 
respect to the grid spacing h. By choosing very high order 
interpolating scaling functions we can get arbitrarily high 
convergence rates. Since convolutions are performed with 
FFT techniques the numerical effort does not increase as 
the order m is increased. The accuracy shown in Figure 3 
for the Martyna and Tuckerman method is the accuracy 
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FIG. 2 Plots of the value of Ko(p) (top) and error of the 
integration defining its value (bottom) of for 1024 integration 
points for all the values of p used in the tensor decomposition 
of 1/r in gaussian functions. 



in the central part of the cube that has 1/8 of the total 
volume of the computational cell. Outside this volume 
errors blow up. So the main disadvantage of this method 
is that a very large computational volume is needed in or- 
der to obtain accurate results in a sufficiently large target 
volume. For this reason the less acurate Hockney method 
is generally prefered in the CPMD program (14). 

A strictly localized charge distribution, i.e. a charge 
distribution that is exactly zero outside a finite volume, 
can not be represented by a finite number of plane waves. 
This is an inherent contradiction in all the plane wave 
methods for the solution of Poisson's equation under free 
boundary conditions. For the test shown in Figure 3 
we used a Gaussian charge distribution whose potential 
can be calculated analytically. The Gaussian was em- 
bedded in a computational cell that was so large that 
the tails of the Gaussian were cut off at an amplitude 
of less than l.e-16. A Gaussian can well be represented 
by a relatively small number of plane waves and so the 
above described problem is not important. For other lo- 
calized charge distributions that are less smooth a finite 
Fourier representation is worse and leads to a spilling of 
the charge density out of the original localization volume. 
This will lead to inaccuracies in the potential. 

Table I shows the required CPU time for a 128 3 prob- 
lem as a function of the number of processors on a Cray 
parallel computer. The parallel version is based on a 




8th 
14th 
20th 
30th 
40th - 
50th 
60th 
100th 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 
Grid step 

FIG. 3 Accuracy comparison between our method with inter- 
polating scaling functions of different orders and the Hockney 
of Martyna-Tuckerman method as implemented in CPMD. 
The accuracy of our method is finally limited by the accuracy 
of the expansion of Eq. 12 with 89 terms. 



parallel 3-dimensional FFT. 



1 


2 


4 


8 


16 


32 


64 


.92 


.55 


.27 


.16 


.11 


.08 


.09 



TABLE I The elapsed time in seconds required on a Cray 
XT3 (based on AMD Opteron processors) to solve Poisson's 
equation on a 128 3 grid as a function of the number of proces- 
sors. Since Poisson's equation is typically solved many times, 
the time for setting up the Kernel is not included. Including 
the set up time of the Kernel increases the total timing by 
about 50 percent, since one additional FFT is needed. 



A package for solving Poisson's equation according 
to the method described here can be downloaded from 
www.unibas.ch/comphys/comphys/SOFTWARE 

In conclusion, we have presented a method that allows 
to obtain the potential of a localized charge distribution 
under free boundary conditions with a 0(N log AT) scal- 
ing in a mathematically clean way. Even though our 
method has the same scaling behaviour as existing plane 
wave methods it is not a plane wave method in the sense 
that neither the charge density nor the potential are ever 
represented by plane waves. Instead interpolating scal- 
ing functions are used for the representation of the charge 
density. 

We acknowledge interesting discussions with Rein- 
hold Schneider and Robert Harrison. This work was 
supported by the European Commission within the 
Sixth Framework Programme through NEST-BigDFT 
(contract no. BigDFT-511815) and by the Swiss Na- 
tional Science Foundation. Research of G.B. was par- 
tially supported by DOE grant DE-FG02-03ER25583, 
DOE/ORNL grant 4000038129 and DARPA/ARO grant 
W911NF-04-1-0281. The timings were performed at the 
CSCS (Swiss Supercomputing Center) in Manno. 
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APPENDIX A 

In the present Appendix we are going to prove Eq. 
(6): Let <p{x) be an interpolating scaling function of 
Deslariers-Dubuc, of the order m, and Pi lt i 2t i 3 be a three- 
dimensional array of constant coefficients. Let, further, 



(x - ii)4>{x - i2)4>(x - i 3 ). (Al) 



Then, 



E isX,™ = J drx h y l *z l *p(r) (A2) 

*i,*2,*3 

if < h, h, h < fn 

This follows from the fact, proven in reference (16) 
that the first m moments of the scaling function obey 
the formula: 

Mi = J cp(x)x l dx = Si, Z = 0,..,m-1 (A3) 
Shift the integration variable, we have 

J <t>{x-j)x l dx = J <t>{t){t+j) l &t = 



p=0 



Then, inserting (Al) into the right side of (A2), we get: 
J dvx h y h z h p{v) = J x h y h z h Pu.Ws x 

11,12,13 

X (j>{x - h)(j>{x - i2)<}>(x - h)&Y = E Ph,i2,i3 X 

Hi»2,»3 

x J x h 4>(x - ii)dx J y h (f)(y - i 2 )dy x 

' 3 <j)(z-i 3 )dz= E Pii^AsA 1 ^ 2 ^ 3 
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